Modification effects of socioeconomic factors on associations between air pollutants and hand, foot, and mouth disease: A multicity time-series study based on heavily polluted areas in the basin area of Sichuan Province, China

Background Hand, foot, and mouth disease (HFMD) is a serious threat among children in China. Some studies have found that air pollution is associated with HFMD incidence, but the results showed heterogeneity. In this study, we aimed to explore the heterogeneity of associations between air pollutants and the number of HFMD cases and to identify significant socioeconomic effect modifiers. Methods We collected daily surveillance data on HFMD cases in those aged less than 15 years, air pollution variables and meteorological variables from 2015 to 2017 in the basin area of Sichuan Province. We also collected socioeconomic indicator data. We conducted a two-stage multicity time-series analysis. In the first stage, we constructed a distributed lag nonlinear model (DLNM) to obtain cumulative exposure-response curves between each air pollutant and the numbers of HFMD cases for every city. In the second stage, we carried out a multivariable meta-regression to merge the estimations in the first stage and to identify significant socioeconomic effect modifiers. Results We found that PM10, NO2 and O3 concentrations were associated with the number of HFMD cases. An inverted V-shaped association between PM10 and the number of HFMD cases was observed. The overall NO2-HFMD association was a hockey-stick shape. For the relationships of PM10, SO2, NO2, O3 and CO with HFMD counts, approximately 58.5%, 48.4%, 51.0%, 55.6% and 52.5% of the heterogeneity could be explained, respectively. The proportion of primary school students, population density, urbanization rate, number of licensed physicians and number of hospital beds explained part of the heterogeneity and modified the relationships. Conclusion Our study explored the heterogeneity of associations between air pollutants and HFMD counts. The proportion of primary school students, population density, urbanization rate, number of licensed physicians and number of hospital beds could modify the relationships. The results can serve as a reference for relevant public health decision making.

was observed. The overall NO 2 -HFMD association was a hockey-stick shape. For the relationships of PM 10 , SO 2 , NO 2 , O 3 and CO with HFMD counts, approximately 58.5%, 48.4%, 51.0%, 55.6% and 52.5% of the heterogeneity could be explained, respectively. The proportion of primary school students, population density, urbanization rate, number of licensed physicians and number of hospital beds explained part of the heterogeneity and modified the relationships.

Introduction
Hand, foot, and mouth disease (HFMD), caused by enteroviruses, is a highly prevalent infectious disease mainly occurring in children [1,2]. In recent decades, it has emerged in many Asian countries [1,[3][4][5][6]. HFMD has a high incidence and imposes a serious disease burden on those in mainland China. From 2008 to 2019, a total of approximately 22.45 million HFMD cases were reported nationwide, ranking first among all notifiable infectious diseases. Although HFMD is usually self-limiting, it may sometimes lead to serious central nervous system and cardiovascular complications and even death [7]. Approximately 75,881 disabilityadjusted life years (DALYs) are caused by HFMD annually [8]. Therefore, HFMD is a severe threat among children. However, only inactivated enterovirus A71 vaccines have been developed, and specific therapies are still lacking [9,10]. Thus, nonpharmaceutical interventions are necessary to reduce the risk of HFMD.
Identifying environmental risk factors is important for the implementation of nonpharmaceutical HFMD interventions. Meteorological factors, such as temperature and relative humidity, have been found to be associated with HFMD in many studies [11][12][13][14]. In addition, considering persistent air pollution [15], some studies have recently suggested that air pollutants are associated with the HFMD incidence [16][17][18][19][20][21][22]. Nonlinear relationships and delayed effects have been shown for most of these associations. However, controversial results have been reported among different studies. For example, a study in Wuhan [16] found that PM 10 exposure increased the risk of HFMD on lag10-11 according to a generalized additive model (GAM). Another study in Ningbo [19] demonstrated that there was no significant PM 10 -HFMD relationship according to a distributed lag nonlinear model (DLNM). A study in Guangxi [23] stated that the PM 10 -HFMD relationship has a positive linearity according to a principal component regression (PCR). Studies in Chengdu [20] and Shenzhen [21] showed that the cumulative risk curves of PM 10 and HFMD were inverted "V"-shaped according to the DLNM. Among these studies, heterogeneity of the correlations and the shapes of the association curves can be observed. This heterogeneity may be attributed to different study methods, the modification effects socioeconomic variables and so on [24]. Exploring heterogeneity can reveal the influence of effect modifiers on air pollutant-HFMD associations, which may guide HFMD prevention strategies. However, all previous studies on air pollutant-HFMD associations were conducted in a single city. Thus, the exploration of heterogeneity is still lacking.
To determine the reasons for the heterogeneity and obtain average air pollutant-HFMD associations, a multiregional study with a common statistical method is needed. Therefore, we conducted a multicity study based on 17 cities in Sichuan Province with a two-stage timeseries analysis. We selected this area because HFMD in Sichuan Province has always been ranked among the top three notifiable infectious diseases. In addition, the selected 17 cities, located in the basin area, have serious air pollution because of rapid industrialization [25,26]. Moreover, different socioeconomic development levels can be found among cities. Two-stage time-series analyses, consisting of a DLNM and meta-analysis, have been conducted in many epidemiological studies to explore modification effects [13,14,27,28]. The DLNM can characterize nonlinear associations between two variables in the lag-response and exposure-response dimensions [29,30]. The meta-analysis can merge city-specific associations and identify effect modifiers of heterogeneity among cities [27,31].
In this work, we carried out a two-stage time-series analysis to explore the heterogeneity of air pollutant-HFMD case associations and identify socioeconomic effect modifiers in the basin area of Sichuan Province. In the first stage, we obtained cumulative exposure-response curves between air pollutant concentrations and the numbers of HFMD cases for every city. In the second stage, we merged the estimations in the first stage and included the city-specific socioeconomic indicators to identify significant effect modifiers. Our results can serve as reference data for the prevention and control of HFMD.

Ethics statement
Our study was approved by the institutional review board of the School of Public Health, Sichuan University. All HFMD surveillance data were collected from Sichuan Center for Disease Control and Prevention. The study methods were carried out in accordance with relevant guidelines and regulations. Our study was constructed at the population level. Therefore, no confidential information was involved in this study and informed consent was not required.

Study area
The 17 selected cities, with a total area of 185,757 km 2 , are located in the basin area of Sichuan Province, China. The study area has a subtropical monsoon climate. The western and southern parts of the basin area are heavily polluted, while the northeastern area is less polluted. The western area has a larger population and more health resources than the northeastern area.
According to previous studies, 4 types of city-specific indicators were included in our study: (1) demographic variables: the population density (people per km 2 ), birth rate (‰) and proportion of primary school students (‰); (2) economic variables: gross domestic product (GDP) per person (CNY), GDP increase (%) and urbanization rate (%); (3) health resource variables: the number of licensed physicians (per 1000 population) and number of hospital beds (per 1000 population); and (4) a traffic variable: the total number of passengers (trips). City-specific socioeconomic data were collected from the China City Statistical Yearbooks from 2015 to 2017, and each indicator for each city was calculated by arithmetic averaging of the indicators for three years.

Statistical analysis
A two-stage multicity time-series analysis was carried out to characterize associations between air pollutant concentrations and the number of HFMD cases and identify significant effect modifiers in the study area from 2015 to 2017. First, to obtain the cumulative exposureresponse curves of the effects of PM 10 , PM 2.5 , SO 2 , NO 2 , CO, and O 3 concentrations on the number of HFMD cases, a DLNM was constructed for each city. Then, to obtain overall cumulative exposure-response curves, the estimations in the first stage were merged by multivariable meta-analysis. Finally, each city-specific socioeconomic indicator was included to identify whether the indicator could explain the heterogeneity between cities and evaluate its modification effect. PM 10 and PM 2.5 concentrations were found to be highly correlated (r s = 0.94) by Spearman rank correlation analysis (S1 Fig), which may lead to multicollinearity. Therefore, considering that both particulate matter and PM 10 were analyzed in several studies, PM 10 , SO 2 , NO 2 , O 3 , and CO were finally included in our study.

First-stage analysis.
To measure the cumulative exposure-response curve of the relationship between each air pollutant and the number of HFMD cases for each city, a DLNM was constructed as follows: where Y t is the number of HFMD cases on day t(t = 1,2,. . .1096); α is the intercept; cb(�) is the cross-basis function; and P it is the i(i = 1,2,. . .5)th pollutant (i = 1,2,. . .5) on day t. Natural cubic splines with 3 degrees of freedom (df) were applied to characterize the lag-response relationship between each air pollutant concentration and the number of HFMD cases, and natural cubic splines with 2 knots were applied to characterize the exposure-response relationship. The values of the knots were the 33 th and 66 th percentiles of each air pollutant concentrations for 17 cities (Knots: PM 10  . Zero to fourteen days was the lag interval for PM 10 , SO 2 , NO 2 , CO, and O 3 . To control for confounding factors, the mean temperature was included by calculating the simple moving weighted averages (SMAs) in the same lag interval as the air pollutants and a natural cubic spline with 3 df. M jt is the j(j = 1,2,. . .4)th other meteorological confounder, which was incorporated by calculating the SMAs in the same lag interval as the air pollutants. A natural cubic spline with 8 df per calendar year was used to control for seasonal and long-term trends. Holiday t is a binary variable indicating whether day t was a holiday. DOW t is an indicator of the day of a week. A t is the autoregressive term to control for autocorrelation. Sensitivity analysis, using the sum of the quasi Akaike Information Criterion (QAIC) of 17 cities as the criterion, was conducted to determine the forms of variables.

Second-stage analysis.
The multivariable meta-regression was constructed as follows: where b y q is a k-dimensional parameter vector representing the estimate of the association between air pollutants and HFMD in the q th city. S q is a k × k dimensional variance-covariance matrix of b y q . U q , a k × kp dimensional block-diagonal matrix, is the Kronecker expansion of p area-specific indicators u q = [u 1 , u 2 . . .u p ] T . When there is no meta-variable, U = I (k) and β denote the overall mean of the parameter vector b y q . S q and Ψ are variance-covariance matrices representing within-group variation and between-group variation, respectively.
The multivariable meta-regression with only an intercept and restricted maximum likelihood estimation (REML) were adopted to pool the city-specific estimates from the first stage to obtain the overall cumulative exposure-response associations between the number of HFMD cases and air pollutant concentrations. The heterogeneity of these associations was measured by the Cochran Q test and I 2 statistic. Then, each city-specific socioeconomic indicator was independently included in the meta-regression to explore whether the indicator could explain the heterogeneity, and the likelihood ratio (LR) test was carried out to identify significant effect modifiers (P�0.05).
All statistical analyses were conducted in R 4.0.3 using the packages dlnm, splines and mvmeta. Table 1 shows the summary characteristics of all the variables in our study. A total of 201,035 HFMD cases in individuals aged less than 15 years in 17 cities were reported from 2015 to 2017. The air pollution concentrations were different among cities (S2 Fig). Regarding PM 10 , PM 2.5 , SO 2 , CO and O 3 , the southern cities had heavier pollution. Regarding NO 2 , Chengdu was the most polluted city. Variations were also found among city-specific socioeconomic characteristics (S3 Table). The eastern cities had a higher proportion of primary school students, and the northeastern cities had smaller populations and fewer health resources. Fig 1 illustrates the seasonal trends of HFMD cases and air pollutants. The number of HFMD cases had two peaks every year, which were concentrated in May-July and November-January. O 3 had a peak concentration in summer and a lower concentration in winter. In contrast, the other air pollutant concentrations and the AQI were higher in winter and lower in summer.

Overall pooled curves of the air pollutant-HFMD relationship
We merged the city-specific estimates for 17 cities to obtain overall pooled cumulative exposure-response curves and 95% confidence intervals (CIs) (Fig 2). Fig 2A indicates that the association between the PM 10 concentration and the number of HFMD cases has an approximately inverted V-shape. The relative risk (RR) increased between 0-66 μg/m 3 and peaked before decreasing. The curve of NO 2 (Fig 2C) was approximately hockey-stick shaped, and the RR reached 4.47 (95% CI: 1.49, 13.40) at an NO 2 concentration of 90 μg/m 3 . For O 3 (Fig 2D), the RR increased between 0-50 μg/m 3 , decreased slightly between 50-108 μg/m 3 , and then increased with increasing O 3 concentrations. For SO 2 ( Fig 2B) and CO (Fig 2E), the 95% CIs of their curves almost contained 1, meaning that the associations were not statistically significant.  City-specific indicators

Heterogeneity and effect modifiers of the air pollutant-HFMD relationship
The pooled results showed that the city-specific cumulative exposure-response curves for each pollutant were heterogeneous. Table 2

Discussion
This study explored the relationships between air pollutant concentrations and the number of HFMD cases and identified significant socioeconomic effect modifiers. Our results revealed that PM 10 , NO 2 and O 3 concentrations were associated with the number of HFMD cases in a heavily polluted area with a high HFMD incidence and different socioeconomic levels. Moreover, we found heterogeneity for each pollutant in the cumulative exposure-response associations between cities. Regarding the relationships of PM 10 , SO 2 , NO 2 , O 3 and CO with the number of HFMD cases, approximately 58.5%, 48.4%, 51.0%, 55.6% and 52.5% of the variations could be explained by differences among different cities, respectively. The proportion of primary school students, population density, urbanization rate, number of licensed physicians and number of hospital beds could partly explain the heterogeneity and modify the relationships. The exact mechanisms of the effects of air pollutant concentrations on HFMD cases remain unclear and need further exploration. Our results could help to enhance the understanding of the relationships between air pollutant concentrations and the number of HFMD cases and provide a reference for relevant public health decision-making.
For the relationships between air pollutant concentrations and the number of HFMD cases, we found that the overall association between PM 10 and the number of HFMD cases had an approximately inverted V-shape. When the PM 10 concentration was 0-66 μg/m 3 , the correlation was positive. One possible explanation is that PM 10 may lead to lung function damage

PLOS NEGLECTED TROPICAL DISEASES
Heterogeneity of air pollutant-hand, foot, and mouth disease associations between cities and an inflammatory response [32]. Another is that viruses may attach to particles in the air, which may enhance the transmission of the HFMD virus [33][34][35]. When the PM 10 concentration exceeded 66 μg/m 3 , we found a negative association, which was consistent with results in Chengdu [20], Shenzhen [21], and Guangxi [23]. This phenomenon may be attributed to early warning measures implemented by the government. "The Sichuan Provincial Heavy Pollution Weather Emergency Plan" [36], published in 2014, provides the following warning: when the daily average AQI exceeds 200 for 3 or more days, kindergartens and primary schools are advised to avoid outdoor activities to effectively reduce children's exposure to air pollution. In addition, people may opt to wear masks, use air purifiers, and close windows on heavy pollution days as a result of increased awareness of the hazards of particulate matter exposure [37,38], which might also reduce the risk of HFMD. We found that a low O 3 concentration had a weakly positive correlation with the number of HFMD cases, which was similar to studies in Guangxi [23] and Ningbo (Gu et al., 2020) and was inconsistent with studies in Guilin (Yu et al., 2019) and Shenzhen (Yan et al., 2019). Some studies have shown that O 3 exposure could compromise epithelial defenses, increase transmucosal permeability [39], induce oxidative damage to cells and the lining fluids of the airways [40], and cause an acute effect on child lung function [41], which may increase the susceptibility of children to HFMD. However, some other studies have indicated that an appropriate O 3 concentration could reduce virus production and stimulate cytokine production, which may reduce the risk of EV71 HFMD infection [42]. This may lead to a nonsignificant association with HFMD when the O 3 concentration exceeds 50 μg/m 3 . Further research is needed to explore the association between O 3 and HFMD cases and the possible mechanism.
The overall relationship between the NO 2 concentration and the number of HFMD cases was approximately hockey-stick shaped, and the RR was 4.47 (95% CI: 1.49, 13.40) at 90 μg/ m 3 , suggesting that high NO 2 concentrations greatly increase the HFMD risk. A study in Hong Kong found that NO 2 may cause peptic ulcer bleeding [43], which might promote the fecal-oral transmission of HFMD. Our results suggest that children should reduce their outdoor exposure on days with high NO 2 concentrations. The overall association between the SO 2 concentration and HFMD cases was nonsignificant, which was inconsistent with positive associations in studies in Ningbo [44], Shenzhen [21], Hefei [45] and Wuhan [16]. Approximately 95% of the SO 2 concentration data in our study were lower than 30 μg/m 3 , which may lead to a wide confidence interval at high SO 2 concentrations. Therefore, our results only suggest that the association between low SO 2 concentrations and HFMD cases was not significant. We are unsure whether the nonsignificant result at high SO 2 concentrations was due to the wide confidence interval or no actual association between SO 2 and HFMD. Further studies are needed in areas with high SO 2 concentrations.
Regarding the heterogeneity of the air pollutant-HFMD associations, we found that demographic compositions and health resources could partly explain the heterogeneity. Among these factors, the proportion of primary school students was an important effect modifier that explained 10.9% and 5.9% of the heterogeneity of the associations of SO 2 and NO 2 concentrations with the number of HFMD cases, respectively. A high proportion of primary school students weakened the positive effect of NO 2 on the number of HFMD cases and made the association between SO 2 and the number of HFMD cases negative. A possible reason is that the AQI and the concentration of most air pollutants, including SO 2 and NO 2 , are usually higher in winter. During this time, families with children are more likely to reduce the time spent on outdoor activities because of cold weather and early warning measures [37]. Therefore, people in regions with large proportions of primary school students may take more protective measures, which may reduce exposure to HFMD among children. However, the number of licensed physicians, the number of hospital beds, the population density and the urbanization rate were significant effect modifiers. Their effect modifications on the associations between pollutants and the number of HFMD cases were slight, and their CIs on the 10th and 90th percentiles of the modifiers were heavily overlapping. Therefore, it is difficult to say whether these sight modification effects have epidemiological implications. However, these results could still provide clues for future research to explore the modification effects on regions with larger differences in these factors. The results also suggest that children in regions with lower health resources particularly need to reduce their exposure to air pollutants.
Our study explored the heterogeneity of the relationships between air pollutant concentrations and the number of HFMD cases and identified significant effect modifiers. We mainly focused on the modification effects of socioeconomic factors, so we selected cities with different levels of economic development and similar climates to avoid interference from meteorological factors. Therefore, a limitation is that meteorological factors may have modification effects, which were not explored in our study. Additional studies in more regions with diverse climates are needed. For the exposure assignment of each city, we chose the nearest neighbor station to obtain the meteorological exposure data and calculated the arithmetic mean of the monitoring stations to obtain the air pollution exposure data. Therefore, another limitation is that misclassification of exposure may occur in some populations. Further research can consider other better approaches, such as land use regression and random forests, to obtain more accurate exposure data. Finally, when we explored the association between one air pollutant and HFMD cases, other air pollutants were also important confounders. However, some air pollutant data were correlated. Although we controlled for collinearity by not including highly correlated variables in the same model, the moderately correlated variables still influenced the results.

Conclusions
Our study found that the air pollutant-HFMD associations between cities showed heterogeneity. Different demographic characteristics and health care resources could partially explain the heterogeneity. Susceptible populations in regions with a smaller proportion of primary school students and fewer medical resources should avoid exposure to SO 2 and NO 2 . Children in regions with lower health resources particularly need to reduce their exposure to air pollutants. The results could be helpful for relevant public health decision-making to prevent HFMD.